Averaging Rotations

Averaging rotations (i.e., quaternions) may not be as simple as you expected. Based on the result of this work, this post goes into the details why a simple weighted average of quaternions does not work and how to implement the real interpolation between quaternions. Part of this post is inspired from Quaternion Weighted Average. Thank Daniel.

The Wrong Way: Simple Weighted Average of Quaternions

Problem Definition

Before I go into the details of today's topic, let me first give a relatively formal definition to the problem I will discuss in the later sections.

(Weighted Average of Rotations) Given a set of rotations (quaternions) and their associated weights where and , we would like to calculate the weighted average rotation such that
(1) is continuous at each with fixed ;
(2) is continuous at each with fixed .

This definition points out for a desired weighted average of rotations, we would like it to maintain robust to slight changes of any of the sub-rotations as well as to slight changes of weights. As we will see soon, however, most popular solutions do not adhere to these two properties, and thus creating the so-called jump rotation phenomenon.

Wrong Way #1: N-lerp as Weighted Average

N-lerp is short for Normalized Interpolation. It's a widely-used, simple quaterion interpolation strategy to blend multiple quaternions. Basically, an N-lerped quaternion can be calculated as:

The norm at the denominator is necessary to make sure is a unit quaternion, thus representing a rotation.

To see what's wrong with the simple N-lerp method, let's look at simple case.

Assume we are interpolating two quaternions, one is a fixed identity quaternion, i.e., . The second is a quaternion simply rotating around the axis, i.e., the world up direction. But in the following examples, the rotation axis of will slightly deviate from the world up direction so as to let you see how the rotation axis and angle will change as rotation changes. We also set .

As a result, the final interpolated quaternion will be:

The pesudo code is given as follows:

FQuat WeightedQuat = Weights[0] * Cameras[0]->GetActorQuat() + Weights[1] * Cameras[1]->GetActorQuat();;
WeightedQuat.Normalize();

// Get rotation axis and angle of the second rotation
FVector Axis;
float Angle;
Cameras[1]->GetActorQuat().ToAxisAndAngle(Axis, Angle);

PrintString("Rotation Axis is: " + Axis + ", Rotation Angle is: " + Angle + ", w is: " + Cameras[1]->GetActorQuat().W);

The right panel shows the rotation of the second quaternion, which represents a camera always looking at the cube from the player's view. I also print 's rotation axis and angle, as well as its component (real part of the quaternion).

We can see that, when the rotation angle is in , the rotation axis is just the world up direction (as I said, with slight deviation) , which can be justified by the z value printed on screen. However, when the rotation angle surpasses , the rotation axis suddenly flips over to the inverse world up direction, as shown by the negative z value. This means that, when the rotation angle is between , the rotation axis will be flipped and the interpolated quaternion will lie in the other half space of the plane.

The following figure shows what is means to say the other half space of the plane.

When rotation angle is in , sweeps through the top right quarter of the plane, and when rotation angle is in , sweeps through the top left quarter of the plane. Merging both, only the top half space of the plane will be covered by . The bottom half space, will never have the chance to be reached. This is where the jump rotation comes from: the jump change of the rotation axis changes the direction the quaternion is interpolated, and thus changes the final interpolated quaternion.

You may ask, why does the rotation axis have to be flipped? Can you just maintain the rotation axis and let the rotation angle span ? Good question! But it will still cause jump rotation at where angle or , in which case the interpolated quaternion will only span the right half space of the plane. The same problem again.

In fact, we can take a closer look at the the we compute . Let's try to derive a closed form of it.

Assume and , where is the rotation angle of quaternion . In the simple case above, , but in fact it can be any unit rotation axis.

Then we have the following derivations:

The final quaternion spans over half of the rotation range of , meaning that it either spans over if , or spans if . This is exactly what we can expect from what we've discussed above.

This can be easily extended to multiple quaternions. More quaternions there are, more prone to slight changes the average quaternion will be.

Wrong Way #2: S-lerp as Weighted Average

If you are familiar with spherical interpolation (S-lerp), you may wonder if it's possible to maintain continuity using S-lerp. Unfortunately, the jump rotation issue will still arise even if you use S-lerp.

Let's use the same example of N-lerp. The first quaternion is , and assume the second quaternion as . The weights are for both and . By the definition of S-lerp (I use the vanilla S-lerp definition for simplicity), we have the S-lerped quaternion :

This result is equivalent to N-lerp, which is not surprising as S-lerp is still manipulating raw quaternions.

We can use the following figure to illustrate why N-lerp and S-lerp are not working in our trivial case:

The angle of the result quaternion is a periodic function with respect to the rotation angle of . Though we stipulate that angle is within , we can actually extend it to by adding a period of . But in the game engine, the value will be reduced into . It's clear that at the point of , the function is discontinuous, i.e., while , and this therefore causes the jump rotation issue.

Why Is N-lerp Not Working

In this section, I will go deeper into the mathematical details of the jump rotation problem, formally proving that simple weighted average with either N-lerp or S-lerp will sometimes (and when) cause the jump rotation problem. This can help us determine in which situations you can safely use N-lerp or S-lerp to average quaternions, which is quite efficient to compute compared to the right method.

Theorem 1: Quaternion Discontinuity of N-lerp

First, we present a lemma which will be used to introduce our theorem.

Lemma 1.1: Given two quaternions and their weights , their N-lerp average quaternion , which is a function of , is discontinuous at rotation angle when , assuming the range angle is where also represents zero rotation.

To prove this lemma, we first rewrite the quaternions as . Then, their weighted average is:

and its norm is:

which gives the explicit expression of :

If is continuous at , or , then it must satisfy that:

with any given and .

We only need to compare to the real part of :

Taking square on both sides and simplifying the equation, we have:

implying that the equality holds only when or or . The solution of is discarded by substituting into the equation above. When or , both sides represent the identity quaternion.

In other words, if , this equation does not hold, and thus, is not continuos at .

We can actually estimate the difference of angle between and . This gives us an intuitive impression of to which extent the rotation will jump. For ease of computation, we take the difference of the squared real parts:

where and are respectively the half rotation angle of quaternion formed by substituting and . Assume , the above equation simplifies to:

This means that the magnitude of rotation jump is related to the cosine of .

Theorem 1: Given quaternions and their where and , the N-lerp average quaternion , which is a function of , s discontinuous at rotation angle when the rotation angle of quaternion is not , assuming the range angle is where also represents zero rotation.

This can be seen by decomposing to:

With Lemma 1.1, Theorem 1 is obvious.

Theorem 1 tells us that the N-lerp style average quaternion will be discontinuous at least at , where quaternion suddenly changes the magnitude of its rotation angle from . For most cases, you cannot make sure the quaternion formed by other sub-quaternions has a rotation angle of zero, so you need to do something to relax Theorem 1 so that there will be less chance for to jump rotate.

The following two figure show us how Theorem 1 performs in practice. Note that in Unreal Engine, the rotation angle range is , so you will see jump rotation at angle rather than we discussed above.

Theorem 2: Weight Discontinuity of N-lerp

When satisfying certain conditions, the average quaternion derived from N-lerp also suffers from discontinuity at particular weight values. We first present the following lemma.

Lemma 2.1: Given two quaternions and their weights , their N-lerp average quaternion , which is function of , is discontinuous at if:
(1) and where ;
(2) and where or .

First we can express as:

and its norm is:

Therefore the explicit form of is:

Considering the real part of :

We then discuss two cases whether or not .

When , we have:

When , it's equivalent to zero and thus is continuous on . To find its discontinuity, we first note that is bounded in as is a quaternion, i.e., a four-dimensional unit vector. Both the denominator and the numerator increase linearly with respect to but at different rates. The only possibility of discontinuity is at some such that , and when traverses from to , the sign of changes.

Let's expand :

which implies:

Through the first and the second equations, we have solutions or . But only the second solution is valid, which applies to all and .

Substituting back into we have:

The proof for the second condition is quite similar. Suppose , we have the following three equations hold:

where . From the first and the second equations, we have two solutions, one is and the other is . Substituting into the third equation, we have three solutions: or or , but none of them give us discontinuity. However, by substituting into the third equation, we have:

Since , the only possible solution is , at which time , and at which time .

Theorem 2: Given quaternions and their weights where and , the N-lerp average quaternion , which is function of , is discontinuous at if:
(1) and where ;
(2) and where or .
where the last quaternions form a new quaternion .

This theorem is obvious when we rewrite as:

The following figure shows that even if sub-quaternions are fixed, changing their weights could also drastically change the final average quaternion, resulting in the jump rotation issue.

The Right Way: Matrix Decomposition for the Largest Eigen-value

The jump rotation problem is rooted in the double cover property of quaternions, i.e., and represent the same rotation, so that quaternions provide a 2:1 mapping of the rotation group. Ideally, changing the sign of should not change the result average, but as we've observed above, it indeed changes the result, which is not what we would want.

We would like a way, serving as a 1:1 mapping to the rotation group, to help us disambiguate from the double cover issue. Matrix, does serve this purpose!

Special Orthogonal Matrix Is 1:1 Mapping to Rotation

First according to Wikipedia, we call the subgroup of orthogonal matrices with determinant +1 the special orthogonal group, denoted by . Then we can easily see that every rotation can be represented uniquely by a matrix . This can be proved by showing that any rotation is a linear transformation of the standard basis , which can be represented by a matrix in . For the inverse, any matrix forms a unique orthonormal basis, which represents a unique rotation upon the standard basis.

Next, we prove that the group of quaternions is 2:1 mapping to . This is obvious because the matrix presentation of a quaternion is:

This is a rotation around vector by a angle of with . It's clear that both and are mapped to .

We can understand quaternions as the process of rotation, whereas conceive matrices as the result of rotation, represented as the basis after rotation.

Problem Reformulation

Right now, the problem becomes: given a set of rotation matrices and their associated weights, how do we compute the average rotation?

A trivial idea is that, now that each matrix represents a basis in , we can interpolate these sets of basis. I will revisit this idea later, but now we can leverage a simpler idea to calculate the average rotation. It's based on a particular definition of average: the average rotation should minimize a weighted sum of the squared Frobenius norms of the rotation matrices, which can be formulated as follows:

Based on this definition, we can derive a very neat closed form of , which maintains a nice property of continuity of interest. I call this method, M-lerp, short for Matrix Interpolation for rotations, and call the final quaternion as .

In this follow subsection, I will go into the details of deriving step by step, basically following this work, but with more elaborations on some of the missing steps.

Solving the Problem

Let's first expand in the form of matrices:

In order to expand , we can rewrite into , assuming any quaternion . Note that is the cross-product matrix:

We can then expand as follows:

From this third equality, we know that .

From the definition of and the expansion of , we can reduce to:

To compute , we have:

which means .

Combining all of these, we have:

Finally, the problem reduces to:

and this corresponds to finding the eigenvector of the largest eigenvalue of .

We can further examine the error quaternion between and , represented by , where is the angle difference between and . On the other hand, can also be interpreted in the form of matrix . Hence,

This implies that the average quaternion that minimizes the weighted sum of Frobenius norms actually minimizes the weighted sum of squared sines of half angles, which is quite intuitive to interpret .

Pseudo Code

There are plenty of algorithms to compute the largest eigenvalue of a symmetric matrix, but the most commonly used one is power iteration, as it's simple to implement and fast in most cases.

I won't go into the details of power iteration. You can easily understand this algorithm with basic linear algebra knowledge.

Giving a set of quaternions and their normalized weights , we can compute the average quaternion as follows:

FQuat AverageQuats(TArray<FQuat>& Quats, TArray<float>& Weights)
{
    FMatrix M = FMatrix(
        0, 0, 0, 0,
        0, 0, 0, 0,
        0, 0, 0, 0,
        0, 0, 0, 0
    );

    for (int i = 0; i < Quats.Num(); ++i)
    {
        FQuat Q = Quats[i];

        M += Weights[i] * FMatrix(
            Q.X * Q.X, Q.X * Q.Y, Q.X * Q.Z, Q.X * Q.W,
            Q.Y * Q.X, Q.Y * Q.Y, Q.Y * Q.Z, Q.Y * Q.W,
            Q.Z * Q.X, Q.Z * Q.Y, Q.Z * Q.Z, Q.Z * Q.W,
            Q.W * Q.X, Q.W * Q.Y, Q.W * Q.Z, Q.W * Q.W
        );
    }

    FVector4 V = PowerIteration(M, FVector(0, 0, 0, 1), 64);
    FQuat Result = FQuat(V.X, V.Y, Y.Z, Y.W);

    return Result;
}

FVector4 PowerIteration(FMatrix& M, FVector4& Init, int Steps, float Eps = 1e-5f)
{
    FVector4 EigenVector = Init;
    float EigenValue = M.Product(EigenVector).X / EigenVector.X;

    for (int i = 0; i < Steps; i++)
    {
        FVector NextVector = M.Product(EigenVector);
        NextVector = NextVector.Normalize();
        NextValue = M.Product(NextVector).X / NextVector.X;

        if (FMath::Abs(EigenValue - NextValue) < Eps)
        {
            break;
        }

        EigenVector = NextVector;
        EigenValue = NextValue;
    }

    return EigenVector;
}

Why Is M-lerp Working

Using the same notion of continuity, we conclude that the matrix decomposition method quite well maintains the property of continuity. This can be formalized as the following theorems.

Theorem 3: Quaternion Continuity of M-lerp

Without loss of generality, we view as a function of , i.e., where . We will show that this function is continuous on .

For ease of illustration, we rewrite as where all variables are fixed except . Then, the average quaternion will be .

We first show that the mapping from a real symmetric matrix to its largest eigenvalue is continuous.

Lemma 3.1: The mapping from a real symmetric matrix to its largest eigenvalue is continuous.

We begin with that fact that the spectral norm of is given by its largest singular value, which is the square root of the largest eigenvalue of . Since is symmetric, we have , and its largest eigenvalue is the squared largest eigenvalue of , and thus the spectral norm of is 's largest eigenvalue.

On the other hand, we know that the spectral norm map is continuous, and thus the map is continuous.

Then, we prove the following lemma.

Lemma 3.2: The mapping from a unit vector to matrix is continuous regarding the Frobenius norm.

Assume is any unit vector and . For and close enough to (in terms of the angle between and ):

As long as , , which proves that the mapping is continuous.

Corollary 3.1: The mapping from a unit vector to matrix is continuous regarding the Frobenius norm, where are fixed coefficients and are fixed unit vectors.

Lemma 3.3: The mapping from a unit vector to the eigenvector of the largest eigenvalue of is continuous. That is, is continuous, where .

In fact, the matrix has only one non-zero eigenvalue, 1, and its eigenvector is . Proof is quite simple:

This means the eigenvector is proportional to . Then, assume and we have:

and any vector parallel to is the eigenvector. For simplicity, we take as the eigenvector. This ends our proof.

Theorem 3: The mapping from a quaternion (4-dimensional unit vector) to the eigenvector of the largest eigenvalue of matrix is continuous.

For ease of elaboration, we rewrite . For any quaternion and close enough quaternion (in terms of their angle ), we know from Lemma 3.1 and Corollary 3.1 that:

where is the largest eigenvalue of and is the largest eigenvalue of . It follows that:

where and are respectively the eigenvector of and , i.e., and .

Since is close enough to , then from Lemma 3.2, there exists a matrix such that and .

On the other hand, we use a rotation matrix to represent the rotation from to , i.e., . Note that is a rotation matrix with being:

It's orthogonal and skew-symmetric: . This matrix is actually the matrix representation of quaternion with where is the angle between and . We will show that as .

Substituting and into the above "converges-to expression":

Rearranging the terms gives us:

According to this expression, we know that as . This implies that , indicating that .

We can compute the Frobenius norm of :

implies , and thus . Proof completes.

Remarks: It seems that this is a special case of a more general theorem. Please refer to Kato, T., Perturbation theory for linear operators. Springer-Verlag, 1976.

The following figure shows how M-lerp works to blend three sub-quaternions. Although you can observe that sometimes the camera rotates fast, it's generally continuous as sub-quaternions change without jump rotation. We can say that the averaged rotation is continuous but not continuous. I do not know if this statement is true. It's only my conjecture.

Theorem 4: Weight Continuity of M-lerp

In this subsection, we will show that the mapping of where is continuous on each . More formally, we have the following theorem:

Theorem 4: The mapping from a real value to the eigenvector of the largest eigenvalue of matrix is continuous.

We rewrite . We can then directly compute the derivative of with respect to :

which is constant over . This implies is continuous on . For any close enough to , and by Lemma 3.1, we know that:

where is the largest eigenvalue of and is the largest eigenvalue of where . It follows that:

where and are respectively the eigenvector of and , i.e., and .

Since is close enough to and by the fact that is continuous, there exists a matrix such that and . You can prove this by calculating the Frobenius norm of .

Then we can follow the same proof process of Theorem 3 to derive the conclusion.

The following figure shows that as the weights change, the averaged quaternion still maintains continuity.

A Simple Example: A Two-Rotation Case

Let's revisit the toy example introduced in the Wrong Way 1 section, i.e., and . Here we change the order of and to align the notation with the one introduced in this section. We further assume so with . The weights of and are respectively and where .

According to our final optimization goal:

we can expclitly obtain as:

We then need to compute the eigenvector of its largest eigenvalue. First, we find its eigenvalues via determinant:

Hence, matrix has three distinct eigenvalues, as long as and : and . The largest eigenvalue, of course, is .

To obtain its corresponding eigenvector, we need to solve the following linear equations:

Solving this equation, we have:

Then , which means is continuous at , and thus on . This shows that the averaged rotation is indeed continuous about both the input rotation and the weight.

The Not-Right Way: Quaternion Flipping

The most common implementation of averaing rotations might be to maintain a running average and flip any incoming sub-rotation if the dot product between and is negative, i.e., it's on the opposite hemisphere to the running average.

FQuat RunningQuat = FQuat::Zero();

for (int i = 0; i < Quats.Num(); ++i)
{
    if (FQuat::Dot(RunningQuat, Quats[i]) >= 0)
    {
        RunningQuat = RunningQuat + Weights[i] * Quats[i];
    }
    else  // Flip
    {
        RunningQuat = RunningQuat - Weights[i] * Quats[i];
    }
}

RunningQuat.Normalize();

In this section, I will show that this method also introduces discontinuity in most cases, but it indeed grants continuity under certain conditions. If these conditions are satisfied in your application, you can safely use this approach to average rotations. We call this method F-lerp.

Theorem 5: Conditional Continuity of F-lerp

We first present a special case of F-lerp to give you an intuitive insight of how F-lerp ensures continuity.

Theorem 5: Given quaternions and their weights with , their L-lerp average rotation is continuous about if for any pair of , their dot product is non-negative, i.e., . In this case, F-lerp reduces to N-lerp.

We explicitly write as:

To ensure its continuity, we must have , which implies:

We only need to compare the real part of both sides:

Simplifying this equation we have:

Notice that when , the corresponding quaternion , and according to the given condition that , we have . Similarly when , we have . This means , and thus . Proof completes.

The following figure shows three rotations that are close enough are averaged without any discontinuity.

A corollary of Theorem 5 is stated as follows:

Corollary 5.1: Given quaternions and their weights with , their L-lerp average rotation in the order of is continuous about for some if and .

Proof is omitted.

Theorem 5 and its Corollary 5.1 tell us that, if either all quaternions are close enough, or only one quaternion is far away from other quaternions that are close enough, F-lerp ensures continuity.

Generalizing Corollary 5.1, we have the following Corollary 5.2:

Corollary 5.2: Given quaternions and their weights with , if can be partitioned into two sets and , where and , and the two sets satisfy (1) , (2) , and (3) , then the result of F-lerp is continuous regarding each quaternion.

Corollary 5.2 tells us that if all quaternions can be divided into two sets where intra-dot product is positive (the quaternions in each set are close enough) and the inter-dot product is negative (quaternions from different sets are far away from each other), F-lerp then ensures continuity.

However, the conditions of Corollary 5.2 are usually hard to achieve in practice, adding to the difficulty of applying F-lerp into real-time scenes. In the following section, I will show that, F-lerp is very prone to the slight changes of some points. These points become the discontinuities.

Theorem 6: Discontinuity of F-lerp

However, in this section, I will show you that F-lerp also introduces discontinuities while eliminating some discontinuities brought by N-lerp. We have the following Theorem 6.

Theorem 6: Given quaternions and their associated weights with , the result of F-lerp has discontinuities which satisfy the equation , assuming and all are fixed and is variable. Here is the quaternion derived from the -th iteration of the F-lerp algorithm.

The intuition behind Theorem 6 is that as is continuous regarding all its elements, there should be some points at which the dot product between and changes its sign in the neighbourhood, i.e., from positive to negative or from negative to positive. This causes a significant difference between the two consecutive quaternions, giving rise to jump rotation. Note that by default, refers to the normalized result of the F-lerp algorithm though normalization is only applied at the final step in the original algorithm.

The key point is to find that makes . According to the definition of and , this gives us:

For example, when , this equation reduces to:

which implies that is a discontinuity regardless of .

We can in fact estimate the difference bwtween the quaternion with which the dot product is positive and the quaternion with which the dot product is negative.

Suppose the discontinuity point is and in the positive neighborhood the dot product is positive, and in the negative neighborhood the dot product is negative. The rotation axis is fixed.

The positive result of F-lerp is (assume for ease of calculation):

Similarly, the negative result of F-lerp is:

It is obvious that taking limits for both and , the results are generally not equivalent:

and

The denominator seems a little crazy, but fortunately we can simplify it. Note that satisfy:

Then we have:

Next we expand the denominator of :

Therefore, we can rewrite and as:

The difference between and is determined by and .

We can use a simple example to illustrate this. Suppose and . This assumption satisfies the equation given at the beginning.

Then, we can easily calculate and as follows:

It is clear that means a rotation around the axis by and is a rotation around the axis by , and the difference is .

The following figure illustrates that discontinuity does not appear at but at a point that makes the dot product zero. Jump rotation occurs.

From Theorem 5, we know that if the dot product is always positive, the result of F-lerp will be continuous. The question, however, arises that is any quaternion that makes the dot product zero always a discontinuity? The answer is yes. From the expressions of and , to have them identical, we must have and , which is impossible for a valid quaternion. This can also be seen from the fact that .

Remarks: I cannot figure out why F-lerp becomes the dominant paradigm to average rotations because -- assume the deduction of Theorem 6 is correct -- there are almost no benefits to flip quaternons based on dot product. If these quaternions are close enough, which means that their dot products are always positive, we can simply employ N-lerp to average them according to Theorem 5. If some of them diverge, there is still high chance to cause discontinuity. The only benefit of F-lerp I can think of is that it indeed sometimes remedies the discontinuity induced by N-lerp. As stated in Theorem 1, could be the discontinuity point, but with F-lerp, this discontinuity point can be eliminated. However, as claimed above, F-lerp may introduce other discontinuities. In such cases, I would prefer M-lerp.

The Not-Wrong Way: Make the Discontinuity Continuous

C-lerp: Circular Mean as a Way of Average

Okay let's now forget about everything and think the very first question: why jump rotation happens? The answer is the illness of angle representation. The angle representation in game engine is or , both representing a full circle. In either case, we cannnot simply use the arithmetic mean to average angles because it's undirectional.

What about transforming the trivial angle representation into some other representation that is continuous regarding the input angle? Trigonometric functions, particularly cosine and sine functions, as we all know, serve this purpose quite well, since and .

Therefore, given a set of angles, we can first map them to the unit circle, i.e., converting from polar coordinates to Cartesian coordinates, compute the arithmetic mean of these points, and them map the averaged point back to the angle representation. This will be the final average angle we want.

More formally, assume are angles, we can use the following formula to calculate the average angle :

If each angle is weighted by , then the average angle would be:

This method is called C-lerp, short for circular interpolation, directly from the concept of Circular Mean.

We can write the pseudo code like:

float SumSinYaw = 0.f, SumCosYaw = 0.f;
float SumSinPitch = 0.f, SumCosPitch = 0.f;
float SumSinRoll = 0.f, SumCosRoll = 0.f;

for (int i = 0; i < Cameras.Num(); ++i)
{
   FRotator Q = Cameras[i]->GetActorRotation();
   
   SumSinYaw   += Weights[i] * FMath::DegSin(Q.Yaw);
   SumCosYaw   += Weights[i] * FMath::DegCos(Q.Yaw);
   SumSinPitch += Weights[i] * FMath::DegSin(Q.Pitch);
   SumCosPitch += Weights[i] * FMath::DegCos(Q.Pitch);
   SumSinRoll  += Weights[i] * FMath::DegSin(Q.Roll);
   SumCosRoll  += Weights[i] * FMath::DegCos(Q.Roll);
}

return FRotator(
   FMath::RadiansToDegrees(FMath::Atan2(SumSinPitch, SumCosPitch)),
   FMath::RadiansToDegrees(FMath::Atan2(SumSinYaw, SumCosYaw)), 
   FMath::RadiansToDegrees(FMath::Atan2(SumSinRoll, SumCosRoll))
);

The following figure gives two examples of how C-lerp works:

First, two angles are mapped to their Cartesian coordiantes and . Then the axis and axis values are independetly averaged as . Last, arctan is applied to convert the average point back to polar coordinates.

Looks perfect, right? However, unfortunately vanilla C-lerp does NOT ensure continuity. As a simple example, take and as a variable. According to C-lerp, the average angle would be:

It is obvious that is discontinuous at . When , whereas as . All the discontinuity is rooted in the discontinuity of the function, though and are continuous.

The following figure shows this discontinuity.

Albeit discontinuous at some points, C-lerp can replace L-lerp and F-lerp in most cases as there is less chance to trigger such discontinuity.

SC-lerp: Smoothed C-lerp

Conclusion

References

Averaging Quaternions
Quaternion Weighted Average
Perturbation Theory for Linear Operators